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Abstract 

A  1/8°  global  version  of  the  Navy  Coastal  Ocean  Model  (NCOM)  is  described  with  details  of  its  formu¬ 
lation,  implementation,  and  configuration  of  the  vertical  coordinate.  NCOM  is  a  baroclinic,  hydrostatic, 
Boussinesq,  free-surface  ocean  model  that  allows  its  vertical  coordinate  to  consist  of  a  coordinates  for 
the  upper  layers  and  z-levels  below  a  user-specified  depth.  This  flexibility  allows  implementation  of  a  hybrid 
cr-z  coordinate  system  that  is  expected  to  mitigate  some  of  the  weaknesses  that  can  be  associated  with  either 
pure  coordinate  option.  For  the  global  NCOM  application,  the  cr-z  coordinate  is  used  to  allow  terrain-fol¬ 
lowing  a  coordinates  in  the  upper  ocean,  providing  better  resolution  and  topographic  fidelity  in  shelf 
regions  where  flow  is  most  sensitive  to  its  representation.  Including  z  coordinates  for  deeper  regions  effi¬ 
ciently  maintains  high  near-surface  vertical  resolution  in  the  open  ocean.  Investigation  into  the  impact 
of  the  selected  coordinate  system  focuses  on  differences  between  atmospherically-forced  free-running  (no 
assimilation)  global  solutions  using  cr-z  and  pure  z  coordinates.  Comparisons  with  independent  tempera¬ 
ture  observations  indicate  that  global  NCOM  using  the  cr-z  coordinate  has  improved  skill  relative  to  its 
z  coordinate  implementation.  Among  other  metrics,  we  show  that  in  comparison  with  time  series  of  surface 
temperature  from  National  Oceanic  Data  Center  (NODC)  buoys,  mostly  located  in  coastal  regions,  root 
mean  squared  differences  (RMSD)  improved  for  63%  and  correlation  improved  for  71%  of  the  stations 
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when  a-z  coordinates  were  used  instead  of  pure  z,  For  the  exclusively  open-ocean  Tropical  Atmosphere- 
Ocean  (TAO)  buoys,  differences  between  the  simulations  were  small,  with  the  a-z  showing  smaller  RMSD 
for  45%  of  the  stations  and  higher  correlation  for  65%  of  the  stations.  Additional  comparisons  using  tem¬ 
perature  profile  observations  further  confirm  a  tendency  for  improved  performance  using  the  hybrid  a-z 
coordinates. 

Published  by  Elsevier  Ltd. 
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1.  Introduction 

A  1/8°  global  version  of  the  Navy  Coastal  Ocean  Model  (NCOM)  has  been  developed  by  the 
Naval  Research  Laboratory  (NRL)  and  transitioned  to  begin  operations  at  the  Naval  Oceano¬ 
graphic  Office  (NAVOCEANO).  NCOM  is  a  result  of  efforts  to  standardize  model  development 
at  NRL,  leading  to  more  efficient  configuration  management  and  reduced  duplication  of  effort  for 
transition  and  support  of  models  implemented  for  Naval  operations.  This  paper  covers  the  formu¬ 
lation  of  NCOM,  its  global  implementation,  and  a  comparison  of  results  arising  from  different 
choices  of  the  vertical  coordinate  system.  Additional  information  on  global  NCOM  is  provided 
by  Rhodes  et  al.  (2002),  Barron  et  al.  (2004),  and  Kara  et  al.  (this  issue). 

Entities  that  are  engaged  in  global  maritime  operations  have  an  interest  in  real-time  mesoscale 
analyses  and  forecasts  over  the  global  ocean.  Properties  of  interest  include  sea  level,  temperature, 
salinity,  and  currents.  From  these,  additional  properties  such  as  transport,  density  and  sound 
speed  may  be  derived.  For  operational  applicability,  this  information  must  be  available  at  short 
notice  anywhere  on  the  globe.  Ideally,  the  data  would  be  provided  with  reasonable  accuracy  and 
meaningful  uncertainty  estimates.  If  higher  resolution  information  is  needed,  the  global  data 
should  support  nesting  of  relocatable  models. 

A  variety  of  resources  are  utilized  to  operationally  provide  global  ocean  information.  Satellite 
observations  are  invaluable  for  analyses  of  the  ocean  surface.  However,  the  relatively  fine  spatial 
and  temporal  resolution  of  satellite  surface  observations  is  not  generally  available  for  the  interior 
of  the  ocean.  Direct  subsurface  measurements  tend  to  be  too  sparse  and  widely  scattered  to  re¬ 
solve  ocean  features  over  most  of  the  ocean.  Ocean  General  Circulation  Models  (OGCMs)  help 
compensate  for  the  sparseness  of  oceanic  observations  by  providing  dynamic  interpolation  of  the 
available  data.  Thus,  they  play  an  important  role  in  representing  the  ocean  component  of  the  cli¬ 
mate  system  on  a  wide  variety  of  temporal  and  spatial  scales. 

There  are  a  variety  of  OGCMs  in  general  use.  One  of  the  most  definitive  ways  to  distinguish 
OGCMs  that  use  finite  differences  in  the  horizontal  is  by  their  choice  of  vertical  discretization. 
Lagrangian  and  isopycnal  layer  models  have  primarily  been  used  in  deep-water  modeling  (Hurl- 
burt  and  Thompson,  1980;  Bleck  et  al.,  1992),  and  in  predicting  climatological  and  interannual  sea 
surface  temperature  with  relatively  low  vertical  resolution  as  well  (Kara  et  al.,  2004,  2003).  Recent 
developments  have  shown  broader  application  and  highlighted  the  benefits  of  isopycnal  layers 
within  hybrid  coordinates  (Bleck,  2002;  Chassignet  et  al.,  2003;  Halliwell,  2004).  These  models 
are  efficient  in  describing  mesoscale  ocean  circulation  in  which  the  mesoscale  dynamics  are  depen¬ 
dent  mainly  upon  a  few  vertical  modes.  The  layer  formulation  of  these  models  helps  reduce 
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spurious  numerical  diffusion  in  the  vertical,  and  isopycnal  models  minimize  spurious  cross-isopyc- 
nal  mixing  caused  by  horizontal  diffusion  when  the  isopycnal  surfaces  are  not  aligned  with  the 
surfaces  defined  by  the  vertical  coordinate.  Layer  models  have  difficulty  handling  overturning  cir¬ 
culations,  and  the  intersection  of  the  layers  with  the  ocean  bottom  or  the  ocean  surface  causes 
complications.  Isopycnal  models  can  encounter  difficulty  in  situations  where  different  water 
masses  mix  and  become  ill-defined,  e.g.,  during  deep-water  convection  or  vertical  mixing  in  shal¬ 
low  water. 

While  a  layer  model  pursues  discretization  according  to  the  properties  and  time-varying  thick¬ 
ness  of  each  layer,  level  models  typically  mark  the  water  at  depths  that  are  fixed  relative  to  the 
bottom,  the  z-level  approach,  or  that  are  a  fixed  fraction  of  the  total  water  depth,  the  a  coordinate 
approach.  Isopycnic  models  are  discrete  in  density,  while  non-adaptive  (z,  a)  coordinate  models 
are  continuous  in  density.  In  general,  a  coordinates  have  been  popular  in  shallow-water  and 
coastal  modeling  because  of  their  ability  to  follow  the  bathymetry  and  thereby  provide  more  con¬ 
sistent  vertical  resolution  within  the  bottom  boundary  layer.  They  also  naturally  provide  in¬ 
creased  vertical  resolution  in  shallow  water.  However,  ocean  models  with  a  coordinates  can 
suffer  significant  errors  in  their  horizontal  pressure  gradient,  advection,  and  diffusion  terms  in  re¬ 
gions  of  abrupt  changes  in  depth  (e.g.,  Martin,  1985;  Deleersnijder  and  Beckers,  1992;  Paul,  1994; 
Ezer  et  al.,  2002). 

Alternatively,  z-level  coordinates  may  be  favored  for  allowing  simple  calculation  of  the  hori¬ 
zontal  pressure  gradient  and  avoiding  truncation  errors  that  can  be  encountered  in  calculating  this 
term  along  steeply  sloping  coordinate  surfaces.  Vertical  resolution  at  a  particular  depth  is  con¬ 
stant,  which  provides  control  of  minimum  vertical  resolution  in  deep-water  regions.  But  z-level 
coordinates  can  suffer  from  inaccuracy  in  representing  the  bottom  depth  if  the  bottom  depth  is 
rounded  to  the  nearest  z-level,  and  they  cannot  provide  locally  increased  resolution  in  the  bottom 
boundary  layer  or  in  shallow  water  without  the  expense  of  increased  resolution  over  the  entire 
domain.  Memory  may  also  be  wasted  accounting  for  inactive  points  below  the  ocean  floor. 

Additional  refinements  to  the  model  formulation  with  either  a  or  z  coordinates  can  mitigate 
their  deficiencies  at  various  levels  of  increased  computational  cost.  A  detailed  examination  of  a 
broad  variety  of  methodologies  which  may  be  employed  is  beyond  the  scope  of  this  paper.  We 
present  results  using  the  hybrid  a-z  coordinate,  an  approach  accommodated  by  NCOM  that  is 
designed  to  capitalize  on  some  of  the  advantages  and  avoid  some  of  the  difficulties  posed  by 
the  traditional  coordinates.  These  results  are  compared  with  NCOM  run  using  z  coordinates  with 
a  free  surface  to  illustrate  some  effects  of  coordinate  system  choices. 

The  purposes  of  NCOM  in  general  and  global  NCOM  in  particular  include  support  of  multiple 
ocean  nests  and  coupling  with  atmospheric  models  to  better  represent  air-sea  interactions,  in 
addition  to  providing  stand-alone  ocean  nowcasts  and  forecasts.  In  its  global  configuration,  which 
is  undergoing  operational  testing  at  NAVOCEANO,  NCOM  will  provide  three-  to  five-day  fore¬ 
casts  and  host  an  embedded  Arctic  ice  model  in  addition  to  supporting  various  higher-resolution 
nested  ocean  models. 

This  paper  is  organized  as  follows:  Section  2  describes  the  basic  physics,  numerics,  and  compu¬ 
tational  procedures  used  in  NCOM.  Section  3  describes  the  model  setup  for  global  simulations. 
Section  4  provides  some  comparisons  of  the  use  of  pure  z-level,  and  a-z  (hybrid)  coordinates 
in  NCOM.  Section  5  discusses  the  impact  of  coordinate  system  choice  on  subsurface  tempera¬ 
tures,  and  Section  6  offers  some  conclusions  and  directions  for  further  research. 
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2.  NCOM  description 

The  physics  and  numerics  of  NCOM  are  based  largely  on  the  Princeton  Ocean  Model  (POM)  as 
described  in  Blumberg  and  Mellor  (1987),  with  some  aspects  from  the  Sigma/Z-level  Model 
(SZM)  (Martin,  2000),  and  with  some  additional  features.  POM  is  a  three-dimensional,  primitive 
equation,  baroclinic,  hydrostatic  and  free  surface  model.  It  uses  an  orthogonal-curvilinear  hori¬ 
zontal  grid,  a  a  (i.e.,  bottom-following)  vertical  grid,  a  split-explicit  treatment  of  the  free  surface, 
Smagorinsky  horizontal  mixing,  and  the  Mellor-Yamada  Level  2.5  turbulence  model  for  vertical 
mixing.  SZM  is  similar  in  many  ways  to  POM,  but  differs  from  POM  in  that  it  uses  a  Cartesian 
horizontal  grid  (the  grid  spacing  in  the  horizontal  is  constant),  a  combined  o/z- level  vertical  grid, 
an  implicit  treatment  of  the  free  surface,  horizontal  eddy  coefficients  calculated  based  on  a  maxi¬ 
mum  grid-cell  Reynolds  number  criteria,  and  vertical  eddy  coefficients  calculated  using  the  Mel¬ 
lor-Yamada  Level  2  turbulence  closure  scheme.  Primitive  equations  and  the  finite  differencing  and 
computational  procedures  used  in  NCOM  are  provided  here,  in  detail. 


2.1.  NCOM  physics 


NCOM  has  a  free  surface  and  is  based  on  the  primitive  equations  and  the  hydrostatic,  Bous- 
sinesq,  and  incompressible  approximations.  The  basic  equations  are  given  in  the  Cartesian  coor¬ 
dinate  system  as  follows: 

Continuity: 


3 u  dv  3vv  _  ^ 

3jc  dy  dz 


(1) 


Momentum: 

I + v  ■  <ro> =-jM+fv+Qu+ VhMm  Vku) + 1  (r- 1)  ■  (2) 

Temperature: 

I + V .  (,n  =  QT + v„kh  v„n + a  g  - + 1  (*„  f ) ,  (4) 

Salinity: 

^  +  V  •  (yS)  =  QS  +  Vh(^hVhS)  + 1  (VH  —  ,  (5) 

Elevation: 

0C_3[(C  +  H)u]  m  +  H]v ] 
dl  dx  dy 


(6) 
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Hydrostatic: 

dp/dz  =  -pg,  (7) 

where 


P  —  P(T,S,z), 
f  —  2iQsin(0), 

;3  ;3  -  3 

V  =  ,&+J^  +  k8i’ 

VH=i|-  +  j^. 

ox  oy 


In  the  equations  above,  Q  is  a  volume  source/sink  term  that  can  be  used  to  specify  source  and 
sink  flows  such  as  river  and  runoff  inflows,  /  is  the  Coriolis  parameter,  Qr  is  the  solar  radiation 
penetrating  the  surface,  and  y  is  the  solar  extinction  profile  (Martin  and  Allard,  1993).  Parameters 
appearing  in  these  equations  and  symbols  throughout  the  text  are  briefly  described  in  Appendix 
A.  The  form  of  these  equations  in  a  coordinates  is  given  in  Blumberg  and  Mellor  (1987,  1983). 
The  form  of  the  a  coordinate  equations  used  in  NCOM  is  similar  to  that  presented  by  Blumberg 
and  Mellor  (1987)  except  that  the  depth  in  the  equations  is  replaced  by  min(//,  zf),  where  zs  is  the 
depth  at  which  the  grid  changes  from  a-  to  z-level  coordinates. 

The  surface  boundary  conditions  for  the  above  equations  are  the  surface  stress  for  the  momen¬ 
tum  equations,  the  surface  heat  flux  for  the  temperature  equation,  and  the  effective  surface  salt 
flux  for  the  salinity  equation.  The  bottom  boundary  conditions  are  the  bottom  drag  for  the 
momentum  equations,  which  is  parameterized  by  a  quadratic  drag  law,  and  zero  flux  for  the  tem¬ 
perature  and  salinity  equations. 

The  horizontal  mixing  coefficients  {Am  and  An)  as  parameterized  in  Section  2.3  can  be  calcu¬ 
lated  with  the  Smagorinsky  (1963)  scheme  or  a  grid-cell  Reynolds  number  scheme  where  the  mix¬ 
ing  coefficients  are  determined  from  a  specified  grid-cell  Reynolds  number.  Minimum  values  for 
the  coefficients  can  be  specified  for  both  schemes. 

The  vertical  mixing  coefficients  KM  and  Kn  are  calculated  using  the  Mellor-Yamada  Level  2 
(Mellor  and  Yamada,  1974)  or  Level  2.5  (Mellor  and  Yamada,  1982)  turbulence  models.  For 
the  Level  2  model,  the  turbulent  length  scale  is  set  to  zero  outside  of  turbulent  layers,  and  within 
turbulent  layers  is  computed  as  (( kzt)~l  +  (&zb)-1)-1,  where  zt  and  zb  are  the  distances  to  the 
upper  and  lower  boundaries  of  the  turbulent  layer,  respectively,  and  k  is  von  Karman’s  constant. 
However,  at  the  surface  and  bottom  the  turbulent  length  scales  for  both  models  asymptote  to  k 
times  the  roughness  scale  at  the  boundary.  The  Level  2  model  can  optionally  include  the  Richard- 
son-number  based  mixing  enhancement  scheme  of  Large  et  al.  (1994),  which  provides  for  weak 
mixing  at  the  edge  of  a  turbulent  boundary  layer  for  Richardson  numbers  above  the  normal  crit¬ 
ical  value  of  about  0.25  up  to  a  value  of  0.7.  For  the  Level  2.5  scheme,  the  surface  flux  of  turbulent 
kinetic  energy  can  be  specified  as  in  Craig  and  Banner  (1994). 

Density  is  computed  using  the  Friedrich  and  Levitus  (1972)  polynomial  approximation  of  the 
equation  of  state  or  the  adaptation  by  Mellor  (1991)  of  the  United  Nations  Educational  Scientific 
and  Cultural  Organization  (UNESCO)  equation  of  state  (Millero  et  al.,  1980;  Millero  and 
Poisson,  1981). 
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2.2.  NCOM  numerics 

The  model  equations  are  solved  on  an  Arakawa  C  grid,  and  the  horizontal  grid  is  orthogonal 
curvilinear  as  in  POM  (Blumberg  and  Herring,  1987).  A  distinguishing  feature  of  NCOM  is  its 
vertical  grid,  which  uses  a  coordinates  from  the  surface  down  to  a  user-specified  depth  (zs)  and 
z-levels  below.  The  a  grid  is  divided  into  layers  as  in  POM,  with  each  cr-layer  being  a  fixed  fraction 
of  the  total  depth  occupied  by  the  cr-layers.  On  the  z-level  grid,  the  bottom  depth  is  rounded  to  the 
nearest  z-level. 

Fig.  1  illustrates  the  different  ways  the  grid  can  be  configured.  One  extreme  is  a  vertical  grid 
with  a  single  cr-layer  at  the  surface  and  z-levels  below  (Fig.  la).  Since  the  model  has  a  free 
surface,  at  least  one  a- layer  is  needed  to  allow  for  changes  in  the  surface  elevation.  If  the 
changes  in  the  surface  elevation  are  large  relative  to  the  grid  resolution  desired  near  the  sur¬ 
face,  a  single  cr-layer  may  not  be  sufficient  to  resolve  the  perturbations  in  the  surface  elevation. 
To  address  this  contingency,  several  cr-layer s  may  be  used  in  the  upper  portion  of  the  water 
column,  allowing  fractional  distribution  of  variations  in  surface  elevation  among  each. 
(Fig.  lb). 

For  water  depths  shallower  than  the  zs  transition  depth  from  cr-layers  to  z-levels,  the  cr-layers 
extend  to  the  bottom  and  become  more  concentrated  as  depth  decreases,  maintaining  spacing  that 
is  a  constant  fraction  of  the  water  column.  Fig.  lc  shows  a  grid  in  which  cr-layers  extend  to  the 
bottom  over  the  shelf,  and  z-levels  represent  the  deeper  water  off  the  shelf.  A  grid  in  which  cr- lay¬ 
ers  are  used  all  the  way  to  the  bottom  everywhere  is  also  possible  (Fig.  Id).  Over  deep  water, 
restricting  cr-layers  to  the  upper  portion  provides  better  control  on  the  uniformity  of  near-surface 
resolution. 


Fig.  1.  Schematic  illustrations  of  vertical  grid  options  available  in  the  Navy  Coastal  Ocean  Model  (NCOM):  (a)  a  single 
cr-layer  at  the  surface  and  z-levels  below,  (b)  several  cr-layers  near  the  surface  with  z-levels  below,  (c)  cr-layers  to  the 
bottom  over  the  shelf  and  above  z-levels  in  deeper  water,  and  (d)  cr-layers  everywhere.  Judicious  partitioning  between 
cr-  and  z-levels  can  mitigate  problems  associated  with  either  pure  coordinate. 
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2.3.  Finite  differencing  in  NCOM 


NCOM  uses  a  staggered  Arakawa  C-grid  with  potential  temperature  (T),  salinity  ( S)  and  den¬ 
sity  p.  The  pressure  (p )  is  defined  at  the  center  of  the  grid  cells,  while  the  velocity  components 
(u,v,w)  are  defined  at  the  center  of  the  faces  of  the  grid  cells.  The  eddy  coefficients 
are  calculated  at  the  staggered  grid  locations;  i.e.,  at  the  same  locations  as  the 

velocities. 

NCOM’s  basic  spatial  averages  and  finite  differences  are  mainly  second  order  with  some  options 
for  higher-order  formulations.  The  spatial  differencing  of  the  horizontal  pressure  gradient  term  on 
the  o  part  of  the  grid  is  evaluated  as  in  Blumberg  and  Mellor  (1987).  There  is  an  option  for  the 
quasi-third  order  upwind  advection  scheme  described  by  Holland  et  al.  (1998)  for  momentum 
and  scalars  and  an  option  for  the  Flux-Corrected  Transport  (FCT)  advection  scheme  (Zalesak, 
1979)  for  scalars,  which  avoids  advective  overshoots.  There  are  also  options  for  fourth-order  eval¬ 
uation  of  the  horizontal  pressure  gradient  (McCalpin,  1994)  and  fourth-order  interpolation  of  the 
Coriolis  terms.  The  higher-order  formulations  are  fairly  efficient  to  use  except  for  the  FCT  advec¬ 
tion  scheme,  which  increases  the  model’s  running  time  by  50%  or  more.  Of  these  higher-order  op¬ 
tions,  only  the  quasi-third  order  upwind  advection  is  employed  within  the  cases  presented  here. 

Temporal  differencing  is  leap-frog  with  an  Asselin  (1972)  filter  to  suppress  timesplitting.  All 
terms  are  treated  explicitly  in  time  except  for  the  solution  for  the  free  surface  and  vertical  diffu¬ 
sion.  In  the  solution  for  the  free  surface,  the  surface  pressure  gradient  terms  in  the  depth-averaged 
momentum  equations  and  the  divergence  terms  in  the  depth-averaged  continuity  equation  are 
evenly  split  between  the  old  and  new  time  levels  to  minimize  the  damping  of  surface  waves.  In 
the  vertical  diffusion  terms,  the  field  being  diffused  is  taken  to  be  fully  at  the  new  time  level  to 
avoid  diffusive  overshoots  and  vertical  gradient  reversals. 

The  free-surface  mode  is  calculated  implicitly;  therefore,  the  surface  pressure  gradients  and  the 
divergence  terms  in  the  surface  elevation  equation  have  a  component  at  the  new  time  level  being 
calculated. 

Horizontal  mixing  is  lagged,  i.e.,  evaluated  at  n  -  1,  which  is  required  for  stability.  Vertical 
mixing  is  fully  implicit,  so  the  vertical  heat  flux  is  calculated  as 


-K 


M-l 

H 


sr"+l 

0Z 


(8) 


Here  the  vertical  eddy  coefficient  is  evaluated  using  the  values  of  the  model  fields  at  n  -  1,  which 
helps  to  avoid  exciting  timesplitting  behavior  in  the  leap-frog  scheme.  Fully  implicit  vertical  mix¬ 
ing  is  needed  to  avoid  spurious  flip-flopping  of  the  vertical  gradients,  which  can  occur  with  a  par¬ 
tially  implicit  scheme  when  the  vertical  eddy  coefficients  become  very  large  in  regions  of  strong 
vertical  mixing.  The  advection,  baroclinic  pressure  gradient  and  Coriolis  terms  are  evaluated  at 
the  central  time  level  («).  The  bottom  drag  is  quadratic,  as  shown  below.  Note  that  to  avoid  time¬ 
splitting,  the  explicit  speed  term  is  calculated  at  the  n  -  1  time  step. 

du 
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2.4.  Computational  procedures 

New  baroclinic  horizontal  velocities  are  calculated  using  an  estimate  of  the  new  surface  eleva¬ 
tion,  and  the  forcing  terms  for  the  baroclinic  velocities  are  vertically  integrated  to  provide  the 
forcing  terms  needed  for  the  depth-averaged  barotropic  equations.  The  depth-averaged  momen¬ 
tum  and  continuity  equations  are  then  implicitly  solved  for  the  new  surface  elevation  and  depth- 
averaged  velocities.  The  new  3D  fields  of  baroclinic  velocities  calculated  in  the  previous  step  are 
corrected  by  adding  a  depth-independent  correction  so  that  their  vertical  mean  agrees  with  the 
new  depth-averaged  velocities.  This  effectively  corrects  the  3D  velocities  for  the  new  surface 
elevation  gradient.  The  velocity  field  that  will  be  used  to  advect  the  scalar  fields  is  calculated 
by  adding  a  depth-independent  correction  to  the  3D  velocity  fields  at  time  level  n,  so  that  the 
depth-average  of  the  advection  velocities  is  consistent  with  the  depth-averaged  continuity  equa¬ 
tion.  New  values  of  the  potential  temperature  ( T  and  S)  are  calculated  using  the  previously  com¬ 
puted  advection  fields. 

The  correction  of  scalar  advection  field  ensures  that  the  velocity  field  used  to  advect  the  scalars 
is  numerically  non-divergent,  which  avoids  spurious  sources  and  sinks  when  using  the  flux  form  of 
numerical  advection.  However,  in  the  flux  formulation  of  the  advection  of  momentum,  momen¬ 
tum  is  not  strictly  conserved  without  an  iterative  approach  or  other  modification.  The  advection 
of  momentum  terms  are  needed  to  calculate  a  new  elevation,  but  a  new  elevation  is  needed  to 
solve  for  the  momentum  advection  in  a  manner  that  is  completely  conservative.  This  is  a  difficulty 
for  all  implicit,  free-surface  models  and  also  for  free-surface  models  that  use  a  separate,  small 
timestep  for  the  free-surface  equations.  Iteration  of  the  solution  of  the  baroclinic  momentum 
and  depth-averaged  equations  to  eliminate  this  inconsistency  is  provided  for  in  the  model,  with 
global  NCOM  configured  for  a  single  iteration.  More  iterations  can  be  used  where  deemed  appro¬ 
priate,  but  in  tests  that  have  been  conducted,  the  effect  of  additional  iterations  to  remove  the  slight 
inconsistency  between  the  momentum  advection  and  the  change  in  surface  elevation  was  not 
significant. 

Another  difficulty  in  the  numerical  calculation  involves  the  partially  implicit  bottom  drag  term. 
If  the  bottom  drag  is  explicit,  the  implicit  vertical  mixing  and  the  bottom  drag  are  numerically 
decoupled  from  the  solution  of  the  depth-averaged  equations.  However,  there  is  no  decoupling 
when  the  bottom  drag  calculation  involves  the  new  velocities.  The  initial,  uncorrected  estimate 
of  the  new  baroclinic  velocities  will  be  involved  in  the  calculation  of  the  bottom  drag,  which  is 
part  of  the  forcing  term  for  the  barotropic  mode.  For  this  reason,  it  is  important  that  the  initial 
calculation  of  the  baroclinic  velocities  (i.e.,  uncorrected  for  the  new  surface  elevation)  be  as  accu¬ 
rate  as  possible. 

A  procedure  that  is  frequently  used  with  a  coordinates,  and  in  the  a  portion  of  NCOM  as  well, 
is  to  subtract  the  mean  horizontal  density  profile  from  the  density  field  before  calculating  the  hor¬ 
izontal  density  gradient  (Blumberg  and  Mellor,  1987),  and  to  subtract  a  spatially  smooth  (e.g., 
climatological  or  horizontally  averaged)  temperature  and  salinity  from  the  temperature  and  salin¬ 
ity  fields  when  calculating  horizontal  diffusion  of  these  fields  (Mellor  and  Blumberg,  1985).  These 
fields  are  calculated  from  the  horizontal  averages  over  the  entire  domain  of  the  background  cli¬ 
matology  and  vary  only  with  depth.  These  steps  are  taken  to  reduce  truncation  error  when  calcu¬ 
lating  the  horizontal  density  gradient,  and  to  reduce  the  vertical  component  of  diffusion  of 
temperature  and  salinity  that  occurs  when  horizontal  diffusion  is  calculated  along  sloping  a 
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surfaces.  These  options  are  used  in  the  global  cases  here  but  may  be  omitted  if  desired  for  a  par¬ 
ticular  application. 


3.  NCOM  setup  for  global  simulations 

The  global  NCOM  simulations  presented  in  this  paper  are  conducted  on  a  2048  x  1280  curvi¬ 
linear  grid,  a  rotated  reprojected  bipolar  grid  with  a  transition  zone  (Murray,  1996)  extending 
from  80°S  to  the  full  Arctic  cap  (Fig.  2).  From  Antarctica  to  «32°N,  the  grid  is  spherical  with 
stretching  to  maintain  a  grid  horizontal  aspect  ratio  near  1  with  resolution  «l/6°  x  1/5°  cos@  (lon¬ 
gitude  x  latitude),  where  6  is  latitude.  This  translates  to  «14  km  grid  spacing,  or  rs1/8°  latitude,  at 
45°S.  South  of  «72°S  the  latitudinal  grid  spacing  remains  constant  to  reduce  the  number  of  and 
computational  effort  allocated  to  grid  cells  around  Antarctica.  North  of  «32°N  the  grid  enters  a 
transition  zone  to  a  bipolar  reprojected  Arctic  cap  with  «47°N  singularities  over  land  in  Canada 
and  Russia,  to  minimize  the  distortion  of  ocean  grid  cells.  Grid  spacing  in  the  Arctic  cap  is 
5-1 5  km. 

The  model  has  a  total  of  40  vertical  material  layers  (41  interfaces):  19  tr-layers  from  the  sur¬ 
face  to  137  m  depth  and  21  z-levels  from  137  m  depth  to  the  maximum  depth  at  5500  m.  This 
configuration  would  be  best  depicted  by  Fig.  lc.  The  vertical  grid  is  logarithmically  stretched 
so  that  the  open  ocean  uppermost  material  level  has  a  rest  depth  thickness  of  1  m.  Under  this 
logarithmic  stretching,  the  interface  at  137  m  is  taken  as  an  approximate  depth  above  or  near 
a  typical  shelf  break  and  thus  appropriate  for  a  a  to  z  transition.  NCOM  uses  realistic  bottom 


Fig.  2.  The  curvilinear  NCOM  grid:  (a)  northern  hemisphere,  and  (b)  southern  hemisphere.  The  transition  zone  from  a 
stretched  spherical  grid  to  a  rotated  reprojected  bipolar  grid  begins  at  32°N. 
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topography  obtained  from  the  2-min  Digital  Bathymetric  Data  Base-Variable  resolution 
(DBDB-V2)  data  (Ko,  2004,  personal  communication).  This  includes  all  the  continental  shelves 
for  depths  >5  m. 

Global  simulations  on  the  grid  size  used  by  global  NCOM  require  a  lot  of  computer  time,  so  an 
emphasis  of  the  NCOM  design  is  to  maintain  a  scalable  and  portable  computer  code  that  runs 
efficiently  on  different  computational  architectures.  A  month-long  simulation  on  the  global  grid 
takes  approximately  7  h  of  wall-clock  time  on  128  processors  of  the  IBM  P3  located  at  NAV- 
OCEANO,  Stennis  Space  Center,  MS. 

Two  global  NCOM  simulations  are  performed  to  examine  some  of  the  implications  of  using 
different  vertical  coordinate  systems:  (1)  a-z  (hybrid)  and  (2)  pure  z-level.  In  a  rest  situation  with 
no  elevation  anomalies,  these  coordinate  systems  match  in  regions  with  bottom  depth  greater  than 
137  m,  where  both  have  aim  thickness  in  the  uppermost  material  level.  Both  simulations  are 
identical  otherwise.  A  pure  a  coordinate  case  is  not  considered  because  our  emphasis  on  the  upper 
ocean  leads  to  a  requirement  for  the  uppermost  layer  thickness  to  be  approximately  1  m  for  pur¬ 
poses  of  potential  coupling  with  global  atmospheric  models.  Under  such  a  restriction,  the  shallow 
water  levels  became  too  thin  in  the  pure  a  case,  at  least  using  the  same  time  step  as  the  other  cases, 
a  requirement  for  an  equitable  comparison. 

It  must  be  emphasized  that  simulations  discussed  in  this  paper  use  date-specific  atmo¬ 
spheric  forcing  only,  i.e.,  none  of  the  simulations  include  any  date-specific  oceanic  data  assim¬ 
ilation  including  SST.  Assimilative  runs  are  discussed  in  a  companion  paper  (Kara  et  al.,  this 
issue).  To  perform  simulations,  the  1/8°  NCOM  is  first  configured  for  the  two  grid  specifica¬ 
tions  {a-z  and  pure  z),  separately.  It  is  then  forced  with  atmospehric  variables  from  the  Fleet 
Numerical  Meteorology  and  Oceanography  Center  (FNMOC)  Navy  Operational  Global 
Atmospheric  Prediction  System  (NOGAPS)  (Rosmond  et  al.,  2002).  Both  NCOM  simulations 
are  run  from  1  January  to  31  March  in  2000,  and  they  start  from  the  same  initial  state.  Be¬ 
fore  performing  these  interannual  simulations,  the  model  was  first  spun-up  for  a  little  over  six 
years  from  a  climatological  initial  state  to  statistical  energy  equilibrium  with  monthly  mean 
wind  stress  and  surface  heat  fluxes  obtained  from  NOGAPS.  The  model  simulations  were 
then  extended  interannually  using  6  hourly  atmospheric  wind  and  thermal  forcing  from 
NOGAPS. 

The  model  forcing  is  calculated  from  the  following  time-varying  atmospheric  fields:  wind 
stress  at  the  sea  surface  (N  m”2),  air  temperature  at  10  m  above  the  sea  surface  (C),  air  mixing 
ratio  at  10  m  above  the  sea  surface  (g  kg-1),  and  net  shortwave  and  net  longwave  radiation  at 
the  sea  surface  (W  m-2).  The  sensible  and  latent  heat  fluxes  (W  m-2)  entering  the  net  surface 
energy  balance  equation  are  strongly  dependent  on  SST  and  are  calculated  every  time  step  using 
the  model  SST  in  bulk  formulations  that  include  effects  of  air-sea  stability  through  the  exchange 
coefficients  (Kara  et  al.,  2002).  The  annual  climatological  SST  cycle  is  built  into  the  model  to  a 
limited  extent,  Including  air  temperature  in  the  formulations  for  latent  and  sensible  heat  flux 
along  with  model  SST  in  the  bulk  formulation  automatically  provides  a  physically  realistic  ten¬ 
dency  towards  the  correct  SST  as  discussed  in  another  OGCM  study  by  Kara  et  al.  (2003). 
Although  radiation  fluxes  also  depend  on  SST  to  some  extent,  these  fluxes  are  obtained  directly 
from  NOGAPS  in  order  to  use  the  atmospheric  cloud  mask.  For  these  experiments,  the  short¬ 
wave  radiation  has  been  filtered  to  remove  the  bias  caused  by  6-h  sampling  but  also  eliminating 
the  diurnal  heating  signal. 
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4.  Evaluation  of  vertical  coordinate  system  choices 

Because  all  the  atmospheric  forcing  is  interannual,  the  SST  predicted  by  NCOM  can  be  com¬ 
pared  with  interannual  SST  obtained  from  buoy  measurements.  Daily  averaged  buoy  SSTs  are 
used  from  three  sources:  (1)  the  Tropical  Atmosphere-Ocean  (TAO)  array  (McPhaden  et  al., 
1998),  (2)  the  Pilot  Research  Moored  Array  (PIRATA)  (Servain  et  al.,  1998),  and  (3)  the  National 
Oceanic  Data  Center  (NODC)  database  (more  details  available  at  www.nodc.noaa.gov/BUOY/ 
buoy.html).  All  the  buoys  report  hourly  SST  measured  at  a  depth  of  m  below  the  sea  surface. 
For  the  model-data  comparisons,  daily-averaged  buoy  SSTs  were  formed.  The  reader  is  referred 
to  Kara  et  al.  (this  issue)  for  a  more  detailed  information  about  the  buoy  data  and  their  use  in  the 
model  validation.  SSTs  from  NCOM  were  also  extracted  at  the  buoy  locations.  Since  the  filtering 
to  minimize  sampling  bias  in  the  thermal  forcing  removes  the  diurnal  forcing  signal  for  these  com¬ 
parisons,  NCOM  in  these  evaluations  does  not  simulate  the  diurnal  cycle  and  there  is  no  need  to 
form  a  daily  average  of  model  SST.  As  mentioned  earlier,  these  are  atmospherically-forced  free- 
running  simulations;  the  model  assimilates  no  data  used  for  model-data  comparisons  nor  any 
other  subsurface  or  SST  observations. 

We  have  chosen  four  buoys  to  demonstrate,  through  simulations  of  daily  SST,  the  impact  of 
using  a-z  versus  pure  z  coordinates  in  the  1/8°  NCOM.  Three  out  of  four  buoys  are  from  NODC, 
while  the  other  is  from  the  TAO  array.  The  NODC  buoys  are  located  near  the  coastal  regions 
(Table  1),  where  possible  differences  between  a-z  versus  pure  z  simulations  of  SST  are  anticipated 
to  arise.  Water  depths  where  these  buoys  are  located  are  also  very  shallow  (<55  m).  In  contrast, 
the  TAO  buoy  in  the  eastern  Equatorial  Ocean  (08°N,  155°W)  is  in  a  location  with  bottom  depth 
5249  m.  This  buoy  is  included  in  the  analyses  to  compare  a-z  and  pure  z  coordinate  systems 
between  coastal  (three  NODC  buoys)  and  open  ocean  (a  TAO  buoy)  locations. 

The  effects  of  running  the  model  with  the  pure  z  coordinate  rather  than  the  a-z  coordinate  are 
clearly  seen  at  the  three  NODC  buoys  (Fig.  3).  The  model  SST  obtained  from  the  pure  z  simula¬ 
tion  tends  to  be  colder  than  that  obtained  from  the  a-z  simulation  on  most  of  the  days  at  (40°N, 
073°W)  and  (42°N,  071  °W)  from  1  January  to  31  March  2000.  In  general,  the  SST  time  series 
from  the  a-z  simulation  follows  buoy  SST  better  than  that  from  the  pure  z  simulations  at  NODC 
buoys,  but  there  is  almost  no  difference  in  SST  obtained  from  the  pure  z  simulation  versus  the  a-z 


Table  1 


Sample  buoy  locations  from  the  National  Oceanic  Data  Center  (NODC)  and  Tropical  Atmosphere  Ocean  (TAO) 
arrays  used  in  the  text 


Buoy  name 

Buoy:  ID 

Buoy  location 

Latitude-longitude  (°) 

Depth  (m) 

Delaware  Bay 

NODC:  44009 

(38°N,  075°W) 

38.45°N-074.70°W 

28 

Long  Island 

NODC:  44025 

(40°N,  073°W) 

40.25°N-073.17°W 

40 

Boston 

NODC:  44013 

(42°N,  071°W) 

42.35°N-070.68°W 

55 

Equator 

TAO:  Pacific 

(08°N,  155°W) 

07.97°N-155.00°W 

5249 

The  list  includes  buoy  name,  buoy  originator  (NODC  or  TAO),  approximate  coordinates,  for  simplicity,  as  used  in  the 
text  and  more  precise  latitude  and  longitude.  Water  depth  where  the  buoy  is  located  is  also  given.  Daily  averages  of  SST 
for  each  buoy  are  constructed  from  hourly  SST  values.  The  buoy  near  the  mouth  of  Delaware  Bay  is  roughly  48  km 
southeast  of  Cape  May,  New  Jersey,  USA,  the  one  near  Long  Island  is  about  61  km  south  of  Islip,  New  York,  USA, 
and  the  northernmost  buoy  is  almost  30  km  east  of  Boston,  Massachusetts,  USA. 
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Fig.  3.  Bottom  topography  formed  from  the  1  min  Digital  Bathymetric  Data  Base-Variable  (DBDB-V)  resolution  data 
set  of  NAVOCEANO.  Two  different  color  bars  are  used  to  emphasize  different  depth  ranges:  deep  water  and  shallow 
water  (<200  m)  (top  left  and  right  panels).  Daily  SST  time  series  obtained  from  buoys  (black),  1/8°  NCOM  using  the 
(7-2  (red),  and  1/8°  NCOM  using  the  pure  z  (cyan)  coordinates  at  four  locations  (three  NODC  buoys  and  one  TAO 
buoy)  over  1  January  2000-31  March  2000.  Also  shown  are  daily  SST  differences  between  buoys  and  1/8°  simulations 
using  the  a-z  and  pure  z  coordinates,  separately.  The  reader  is  referred  to  the  text  for  detailed  buoy  specifics.  Note  that 
there  is  no  data  assimilation  in  the  NCOM  simulations. 
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simulation  at  the  TAO  buoy  located  over  the  open  ocean.  As  o-z  vertical  coordinates  designate 
the  shallow  depths  to  be  represented  by  the  a  coordinates,  the  between  the  pure  z  and  o-z  are  usu¬ 
ally  expected  near  coastal  boundaries  at  shallow  depths. 

We  further  analyze  the  model  results  to  demonstrate  the  impact  of  vertical  coordinate  choices 
in  the  global  NCOM  simulations.  Several  statistical  metrics  are  used  for  comparing  the  SST  time 
series  between  buoy  and  1/8°  NCOM.  Let  X,{i  =  1,2, ...,«)  be  the  set  of  n  buoy  (i.e.,  reference) 
values,  and  let  Y,{i  —  1,2, . . .,  n)  be  the  set  of  corresponding  NCOM  estimates.  Also  let  X{  Y)  and 
oxioY)  be  the  mean  and  standard  deviations  of  the  reference  (estimate)  values,  respectively.  Here, 
n  is  91  days,  spanning  the  1/8°  NCOM  simulation  performed  from  1  January  2000  to  31  March 
2000.  Following  Murphy  (1988),  the  statistical  metrics  used  for  model-data  comparisons  are  as 
follows: 

ME  =Y-X,  (11) 


RMSD  = 


'  i=i 


l'/2 


1  ^ 


R^-J^iXi-XW-Y)/^), 

n  i=i 


(12) 

(13) 


SS  =  1  -  RMSD2/4,  (14) 

where  ME  is  the  bias  or  annual  mean  difference,  RMSD  is  the  root-mean-square  difference,  R  is 
the  correlation  coefficient,  and  SS  is  the  skill  score  (see  also  Kara  et  al.,  this  issue). 

Table  2  gives  statistical  model-data  comparisons  between  the  1/8°  NCOM  and  buoy  SST  time 
series  at  the  four  locations  shown  in  Fig.  3.  ME  values  (i.e.,  annual  mean  SST  biases)  are  generally 


Table  2 


Statistical  verification  of  daily  SST  time  series  between  buoy  and  1/8°  NCOM  at  four  locations 


Buoy 

1/8°  NCOM 

RMSD  (°C) 

ME  (°C) 

0BUOY  (°C) 

^NCOM  (°C) 

R 

SS 

(38°N,  075°W) 

Hybrid  o-z 

1.06 

-0.03 

1.73 

2.49 

0.94 

0.63 

Pure  z 

1.29 

0.23 

1.73 

2.66 

0.92 

0.44 

(40°N,  073°W) 

Hybrid  cr-z 

1.15 

-0.81 

1.69 

1.95 

0.91 

0.54 

Pure  z 

1.55 

-1.34 

1.69 

2.15 

0.95 

0.16 

(42°N,  071°W) 

Hybrid  cr-z 

0.58 

-0.21 

1.11 

1.35 

0.92 

0.72 

Pure  z 

0.88 

-0.61 

1.11 

1.54 

0.94 

0.36 

(08°N,  155°W) 

Hybrid  o-z 

0.24 

-0.11 

0.44 

0.52 

0.91 

Pure  z 

0.26 

-0.11 

0.44 

0.53 

0.90 

0.66 

The  number  of  days  used  in  the  statistical  analysis  is  91  (from  1  January  2000  to  31  March  2000).  Statistics  include  root- 
mean-square  difference  (RMSD),  mean  error  (ME),  standard  deviation  for  buoy  (<tBuoy).  standard  deviation  for 
1/8°  NCOM  (<tNcom)>  correlation  coefficient  (JR),  and  skill  score  (SS).  SS  values  of  1  indicate  perfect  simulations.  Note 
that  all  global  NCOM  simulations  (o-z  and  pure  z)  herein  are  atmospherically-forced  free-running  and  therefore 
include  no  assimilation  of  any  SST. 
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small  (close  to  zero)  between  the  model  predicted  SSTs  and  buoy  SSTs  at  all  locations.  This  is  true 
whether  1/8°  NCOM  simulation  is  performed  with  the  a-z  coordinate  or  pure  z  coordinate.  While 
there  is  no  clear  systematic  bias  at  all  buoy  locations,  at  these  locations  during  three  months  gen¬ 
erally  favoring  cooling,  the  model  in  either  configuration  tends  to  have  a  cold  bias.  The  negative 
ME  values  indicate  that  1/8°  NCOM  slightly  underestimates  SST,  i.e.,  the  model  SST  is  slightly 
cooler  than  the  buoy  SST.  Examination  of  RMS  SST  difference  values  clearly  reveals  that  1/8° 
NCOM  with  the  a-z  coordinate  is  superior  to  the  pure  z  coordinate  case  for  predictions  of 
SST,  although  the  deep  water  station  at  (8°N,  155  °W)  shows  almost  no  difference.  The  RMSD 
with  respect  to  buoy  increases  by  35%  (from  1.15  to  1.55  °C)  at  (40°N,  073  °W),  and  by  52%  (from 
0.58  to  0.88  °C)  at  (42°N,  071°W)  when  using  the  pure  z  instead  of  a-z  coordinate  in  the  model. 
While  it  is  interesting  to  note  at  these  two  locations  that  the  R  values  are  smaller  for  the  a-z  sim¬ 
ulation  in  comparison  to  the  pure  z  simulation,  the  differences  are  not  statistically  significant  from 
each  other  at  the  95%  confidence  interval.  Essentially,  the  model  is  able  to  the  capture  the  phase  of 
SST  variability  quite  well  because  the  R  values  are  high  (>0.90)  at  all  locations  and  cases.  Thus, 
for  the  cases  examined  here,  the  use  of  a-z  instead  of  pure  z  coordinates  has  significant  influence 
on  the  bias,  both  of  the  mean  and  of  the  variance,  but  has  limited  effect  on  the  correlation. 

In  addition  to  RMSD,  the  non-dimensional  SS  is  also  used  for  model-data  comparisons.  The 
reason  for  considering  SS  values  in  examining  model  performance  is  that  SST  biases  are  taken 
into  account  in  the  RMS  differences,  but  the  latter  can  be  small  where  SS  and  R  are  poor  because 
of  a  small  amplitude  of  seasonal  cycle  (or  SST  variations  from  day  to  day)  at  some  locations.  In 
addition,  SST  standard  deviation  is  not  the  same  at  each  buoy,  so  a  non-dimensional  metric  (i.e., 
SS  as  used  here)  is  useful  for  making  fair  comparison  of  performance  among  different  buoys.  It  is 
noted  that  the  model  must  have  SS  >0  to  have  skill,  and  SS  =  1  is  perfect  skill  (see  also  Kara  et  al., 
this  issue).  The  model  success  in  predicting  daily  SST  is  also  evidenced  by  SS  values  which  are 
positive.  Comparing  relative  performance  using  the  different  coordinate  systems,  a  substantial 
improvement  is  noted  in  the  model  simulations  performed  with  the  a-z  coordinate  over  those  per¬ 
formed  with  the  pure  z  coordinate.  For  example,  the  SS  value  at  (40°N,  073°W)  increases  more 
than  a  factor  of  three  (from  0.16  to  0.54  °C)  when  1/8°  NCOM  was  run  with  the  a-z  versus  pure 
z  coordinate.  The  improvement  in  SS  values  is  evident  at  all  of  these  locations  when  using  the  a-z 
coordinate  in  model  simulation.  The  a-z  coordinate  case  produces  SS  above  0.5  at  all  four  buoy 
locations,  while  SS  in  the  z  coordinate  case  exceeds  0.5  only  once. 

Overall  performance  of  1/8°  NCOM  simulations  using  either  a-z  or  pure  z  coordinate  is  exam¬ 
ined  at  all  available  TAO,  NODC  and  PIRATA  buoy  locations  (Fig.  4).  The  variety  of  the  buoy 
locations  improves  our  ability  to  discern  impact  of  coordinate  system  choice  and  determine  which 
provides  superior  performance  in  this  application.  All  TAO  and  PIRATA  buoys  are  fairly  distant 
from  coastal  regions,  mostly  equatorial  Pacific  and  Atlantic,  as  seen  from  Appendix  B.  Results 
using  either  a-z  or  pure  z  simulations  are  expected  to  be  similar  at  these  deep-water,  open  ocean 
stations.  On  the  other  hand,  most  NODC  buoys  are  near  the  coastal  boundaries  (Appendix  C) 
where  the  a-z  coordinate  is  expected  to  provide  superior  performance.  There  are  a  total  of  81 
buoys  during  the  three-month  period:  40  TAO,  38  NODC  and  3  PIRATA  buoys,  located  in  dif¬ 
ferent  regions  of  the  global  ocean.  A  few  data  voids  existing  in  buoy  SST  time  series  are  filled 
using  a  cubic  interpolation  in  time. 

Daily  SST  time  series  from  each  buoy  are  compared  with  series  obtained  from  the  model  sim¬ 
ulations.  Statistical  metrics  are  calculated  over  the  3-month  period  from  1  January  2000  to  31 
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Fig.  4.  Locations  of  NODC,  PIRATA  and  TAO  used  for  model-data  comparisons  in  2000. 

March  2000,  similarly  to  the  analyses  presented  earlier  (see  Table  2).  In  particular,  the  RMSD 
value  is  calculated  for  SST  time  series  between  the  buoy  and  1/8°  NCOM  using  the  a-z  coordinate 
at  each  NODC  buoy  location,  and  the  same  process  is  repeated  or  SST  time  series  between  buoy 
and  1/8°  NCOM  using  the  pure  2  coordinate  (Fig.  5).  There  are  clearly  large  differences  in  the 
RMSD  at  some  of  the  buoy  locations,  depending  the  coordinate  system  used  in  the  model  simu¬ 
lation.  To  better  analyze  RMS  SST  difference  with  respect  to  daily  SST  measurements,  NODC 


Fig.  5.  SST  root-mean-square  difference  (RMSD)  between  buoy  and  1/8°  NCOM  at  NODC  buoys  over  1  January 
200(1-31  March  2000.  Values  are  calculated  for  0-2  and  pure  2  simulations  based  on  the  daily  SST  time  series  during  the 
3-month  period.  In  the  x-axis,  D  denotes  a  buoy  located  in  relatively  deep  water,  while  S  denotes  a  buoy  in 
comparatively  shallow  regions.  The  model  is  atmospherically-forced  free-running  and  has  no  assimilation  of  any  SST  or 
other  data.  The  reader  is  referred  to  the  text  for  details  of  statistical  metrics  and  model  simulations. 
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buoys  are  separated  into  two  categories  as  shown  in  the  figure:  (1)  ones  where  the  depths  of  water 
are  greater  than  400  m  (denoted  as  D),  and  (2)  ones  where  the  depths  of  water  are  equal  shallower 
than  400  m  (denoted  as  S).  While  there  are  almost  no  differences  in  RMS  SST  differences  for 
buoys  located  over  deep  regions,  the  same  is  not  true  for  the  ones  near  the  coastal  regions. 

To  highlight  relative  performance  in  modeling  SST,  we  subtract  the  RMSD  for  the  z  case  from 
the  RMSD  for  the  a-z  ,  where  the  individual  RMSD  are  relative  to  the  buoy  observations  as  in 
Fig.  5.  Thus  a  positive  difference  indicates  that  the  z  configuration  produced  closer  agreement  with 
observations  than  did  the  a-z  version.  The  results  in  Fig.  6  tend  negative,  indicating  superior  per¬ 
formance  using  the  a-z  grid.  In  fact,  the  RMS  SST  difference  can  increase  as  much  as  «0.5  °C  at 
(29°N,  080°W)  and  (33°N,  079°W)  when  using  the  pure  z  coordinate  as  opposed  to  the  a-z  coor¬ 
dinate.  While  these  differences  may  appear  to  be  small,  it  should  be  noted  that  they  correspond  to 
a  RMSD  increase  of  44%  (from  1.06  to  1.53  °C)  at  the  former  and  25%  (from  1.98  to  2.47  °C)  at 
the  latter  location.  RMSD  with  respect  to  observations  is  relatively  large  even  for  the  a-z  simu¬ 
lation  at  (33°N,  079°W).  A  variety  of  factors  beyond  the  vertical  coordinate  may  be  impacting  the 
model  performance,  and  the  impact  of  these  may  be  relatively  insensitive  to  coordinate  choice. 
NCOM  suffers  from  atmospheric  forcing  that  contains  errors  due  to  land  contamination  near 
the  coastal  regions,  as  in  other  OGCM  studies  (e.g.,  Kara  et  al.,  2005);  this  would  primarily  be 
reflected  in  comparisons  with  NODC  buoys.  The  atmospheric  fields  used  in  the  model  simula¬ 
tions,  NOGAPS  in  this  case,  is  on  a  1°  grid,  but  NCOM  is  on  a  much  finer  grid  with  resolution 
near  1/8°.  Interpolation  of  atmospheric  forcing  fields  from  NOGAPS  to  the  model  grid  usually 
results  in  errors,  especially  near  the  coastal  boundaries,  as  in  (33°N,  079°W)  because  the  1°  grid 
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Fig.  6.  The  difference  in  RMSD  values  shown  in  Fig.  5.  It  should  be  noted  that  the  RMSD  with  respect  to  the  buoy  is 
first  obtained  for  each  simulation,  separately,  and  then  subtracted  from  each  other.  In  other  words,  in  the  figure  legend, 
RMSD  [Hybrid(<r-z)-pure  z]  represents  the  RMSD  with  respect  to  the  buoy  for  the  pure  z  simulation  subtracted  from 
the  for  the  a-z  simulation.  There  are  a  total  of  24  out  of  38  NODC  buoys  where  the  the  RMS  SST  difference  values  for 
the  pure  z  simulation  are  larger  than  those  for  the  a-z  simulation. 
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from  NOGAPS  cannot  distinguish  a  precise  transition  between  land  and  sea,  significant  for  heat 
fluxes,  for  example.  Further  details  of  this  problem  are  beyond  the  scope  of  this  paper. 

Unlike  the  RMS  SST  difference,  the  R  values  between  the  NODC  buoy  and  1/8°  NCOM  SST 
time  series  are  generally  similar  to  each  other  for  a-z  and  pure  z  simulations  (Fig.  7)  with  the 
exception  of  a  few  locations.  Overall,  the  R  values  obtained  from  the  a-z  simulation  are  higher 
than  those  from  the  pure  z  simulation.  There  are  only  1 1  out  of  38  NODC  buoys  for  which 
the  R  values  from  the  a-z  simulation  («l/4  of  the  buoys)  are  lower  than  those  from  the  pure  z 
simulation  (Fig.  8).  The  maximum  difference  in  R,  0.14,  is  noted  at  37°N,  122°W.  While  the  R 
is  0.55  for  the  pure  z  simulation,  it  improves  to  0.69  for  the  a-z  simulation.  The  latter  R  value 
is  statistically  significant  in  comparison  to  a  R  value  of  0.6,  while  the  former  is  not. 

Finally,  we  examine  how  the  model  simulations  with  a-z  and  pure  z  coordinates  perform  for 
predictions  of  SST  over  the  open  ocean.  Such  an  investigation  can  be  made  using  TAO  and  PIR- 
ATA  buoys  located  in  the  equatorial  Pacific  and  Atlantic  Ocean,  where  the  depth  of  water  is 
>4000  m  at  most  observation  sites.  Unlike  the  NODC  buoys,  the  differences  in  RMS  SST  differ¬ 
ence  for  these  locations  do  not  indicate  that  the  simulation  using  the  a-z  is  superior,  in  terms  of 
SST  simulation,  to  the  case  using  the  pure  z  coordinate  (Fig.  9).  There  are  18  out  of  40  TAO  buoys 
where  the  RMS  SST  difference  with  respect  to  the  buoy  from  the  a-z  simulation  is  lower  than  the 
one  from  the  pure  z  simulation.  Thus,  1/8°  NCOM  performs  relatively  well  at  approximately  half 
of  the  buoy  locations  when  using  the  a-z  simulation,  while  use  of  the  pure  z  coordinate  has  better 
success  at  the  other  buoy  locations. 

Essentially,  one  must  note  that  the  differences  in  RMS  SST  difference  at  the  TAO  locations  are 
much  smaller  than  those  at  the  NODC  locations  (see  Fig.  6).  They  are  generally  close  to  zero  and 
always  within  the  small  range  of  -0.2  to  0.2  °C.  Such  RMS  SST  difference  values  confirm  the  fact 


Fig.  7.  The  correlation  coefficients  (R)  are  between  the  NODC  buoy  observations  and  1/8°  NCOM  SST  time  series 
over  1  January  2000-31  March  2000,  shown  separately  for  the  o-z  and  pure  z  simulations.  Below  the  x-axis,  D  denotes 
a  buoy  located  in  deep  water,  and  S  denotes  a  buoy  near  the  shallow  regions. 
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Fig.  8.  The  difference  in  R  values  shown  in  Fig.  7.  In  the  figure  legend,  R  [Hybrid(<r-z)-pure  z]  represents  the  R  with 
respect  to  the  buoy  for  the  pure  z  simulation  subtracted  from  the  R  for  the  o-z  simulation.  There  are  a  total  of  1 1  out  of 
38  NODC  buoys  where  the  the  SST  correlation  coefficient  for  the  pure  z  simulation  is  higher  than  the  correlation 
coefficient  for  the  <r-z  simulation. 


Fig.  9.  The  difference  in  RMSD  values  at  TAO  buoys.  In  the  figure  legend,  RMSD  [Hybrid(<r-z)-pure  z]  represents  the 
RMSD  with  respect  to  the  buoy  for  the  pure  z  simulation  subtracted  from  the  for  the  o-z  simulation.  There  are  a  total 
of  22  out  of  40  TAO  buoys  where  the  the  RMS  SST  difference  values  for  the  pure  z  simulation  are  lower  than  those  for 
the  <t-z  simulation. 


that  in  the  open  ocean,  there  is  not  a  clear  advantage  of  the  a-z  coordinate  over  the  pure  r  coor¬ 
dinate.  Similar  to  the  differences  in  RMS  SST  difference,  we  note  that  the  R  values  do  not  differ 
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Fig.  10.  The  difference  in  R  values  at  TAO  buoys.  In  the  figure  legend,  R  [Hybrid(<r-z)-pure  z]  represents  the  R  with 
respect  to  the  buoy  for  the  pure  z  simulation  subtracted  from  the  R  for  the  o-z  simulation.  There  are  a  total  of  14  out  of 
40  TAO  buoys  where  the  the  SST  correlation  coefficient  for  the  pure  z  simulation  are  higher  than  those  for  the  a-z 
simulation. 


significantly  from  each  other  when  using  either  coordinate  configuration  at  the  TAO  buoy  loca¬ 
tions  (Fig.  10).  The  R  values  with  respect  to  buoy  from  the  a-z  simulation  are  higher  than  those 
from  the  pure  z  simulation  at  65%  (26  out  of  40)  of  the  TAO  buoy  locations  as  implied  by  the 
differences  in  SST  correlation  coefficient.  On  the  other  hand,  some  of  the  R  values  from  both  sim¬ 
ulations  are  not  statistically  different  from  a  correlation  coefficient  of  0.6,  and  differences  in  most 
of  the  R  values  are  usually  very  small  (<0.05),  resulting  in  a  conclusion  that  the  coordinate  differ¬ 
ences  do  not  have  a  significant  effect  in  deep  water  on  the  shape  of  the  daily  SST  time  series.  The 
magnitude  of  the  R  value  differences  exceeds  0.1  only  at  four  TAO  buoys.  The  model  SST  results 
were  also  analyzed  at  3  PIRATA  buoys,  (0°N,  023°W),  (0°N,  035°W)  and  (15°N,  038°W),  where 
the  SST  time  series  were  available  over  1  January  2000-31  March  2000.  There  were  almost  no  dif¬ 
ferences  between  the  statistical  results  for  these  locations  using  the  a-z  coordinate  as  opposed  to 
the  pure  z  coordinate  in  the  model  simulations  (not  shown),  fn  fact,  the  SST  time  series  were  iden¬ 
tical  from  both  simulations  at  (0°N,  035°W). 

5.  Impact  of  coordinate  system  choice  on  subsurface  temperatures 

Following  the  discussions  presented  in  the  earlier  section,  we  now  examine  how  the  choice  of 
vertical  coordinate  system  (a-z  or  pure  z)  affects  subsurface  temperatures.  Fig.  1 1  shows  cross-sec¬ 
tions  of  subsurface  temperatures  obtained  from  model  simulations  performed  with  a-z  and  pure  z 
coordinates  at  40.25°NN.  This  latitude  for  the  cross-section  analysis  is  chosen  because  there  are 
large  differences  in  SS  values  at  (40°N,  073°W),  as  presented  earlier. 

Differences  are  evident  in  a  comparison  of  the  mean  sections.  One  stark  difference  is  in  the  devel¬ 
opment  of  cold  core  along  the  coast,  which  is  evident  in  the  a-z  simulation  but  much  stronger,  on 
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Fig.  11.  Comparisons  of  vertical  temperature  sections  off  the  eastern  coast  of  the  United  States  at  40.25°N.  The  o-z  and 
pure  z  simulations  show  significant  differences,  with  the  o-z  showing  better  agreement  with  a  time-series  of  buoy  SST 
observation  (Fig.  3).  For  plotting  purposes,  a  dotted  land  mask  is  superimposed  on  filled  temperature  fields;  it  does  not 
reflect  the  discrete  steps  of  the  z-level  bathymetry  that  contribute  to  some  of  the  differences,  (a)  o-z  simulation:  January- 
March  2000  mean,  (b)  o-z  simulation  on  1 1  February  2000,  (c)  o-z  simulation  on  3  March  2000,  (d)  pure  z  simulation: 
January-March  2000  mean,  (e)  pure  z  simulation  on  1 1  February  2000,  (f)  pure  z  simulation  on  3  March  2000. 


the  order  of  1  °C  colder,  in  the  z  case.  Comparison  of  SST  at  (40°N,  073°W)  in  Fig.  3  shows  a  cold 
bias  in  the  z  simulation,  suggesting  that  the  cold  core  in  z  is  too  strong  and  better  represented  by 
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the  a-z  configuration.  Higher  vertical  resolution  in  the  hybrid  coordinate  version  may  play  a  role 
in  its  relatively  better  performance;  for  the  depths  shallower  than  40  m  in  this  zone,  the  a-z  grid 
has  19  levels  while  the  pure  z  version  is  limited  to  no  more  than  12  levels.  The  z  grid  also  shows 
some  evidence  of  warm  patches  along  the  bottom  that  may  be  a  result  of  the  truncated  bathym¬ 
etry.  The  z  case  mean  contains  stronger  horizontal  gradients  of  temperature  inshore,  reflecting  the 
stronger  cold  filament  along  the  coast,  while  the  a-z  simulation  has  stronger  horizontal  gradients 
near  the  100  m  isobath. 

Because  subsurface  temperature  time  series  are  not  available  from  NODC  buoys,  where  the 
largest  SST  differences  are  seen  between  the  a-z  or  pure  z  simulations,  SST  is  used  to  guide  expec¬ 
tations  regarding  the  accuracy  of  subsurface  temperatures  in  the  cross-section  analysis.  Two  spe¬ 
cific  days  are  considered:  1 1  February  2000  and  3  March  2000.  As  indicated  by  Fig.  3,  the  largest 
SST  difference  during  the  3-month  period  at  (40°N,  073°W)  occurs  on  1 1  February  2000.  On  this 
day,  the  SST  from  the  a-z  simulation  agrees  closely  with  the  observation,  while  the  z  case  shows  a 
significant  cold  bias.  We  take  this  to  imply  that  subsurface  temperatures  along  this  section  on  11 
February  2000  are  likely  better  represented  by  the  a-z  simulation  than  by  the  pure  z.  This  com¬ 
parison  is  additional  evidence  that  the  z  configuration  overestimates  the  size  and  strength  of  the 
cold  core  along  the  coast.  At  the  other  extreme,  the  largest  errors  during  the  3-month  period  at 
this  location  occur  in  early  to  mid  March.  On  3  March  2000  the  models  are  in  fairly  close  agree¬ 
ment  with  each  other  but  are  significantly  colder  than  the  observation,  by  ^1 .8  °C  for  the  a-z  sim¬ 
ulation  and  «2.3  °C  for  the  pure  z  case.  Since  the  model  SST  difference  is  comparatively  small  on 
this  date,  we  anticipate  that  subsurface  temperatures  at  (40°N,  073°W)  should  be  similar,  at  least 
up  to  some  extent  near  the  sea  surface.  This  is  the  case  inshore,  but  differences  are  present  at  depth 
along  the  100  m  isobath.  The  a-z  simulation  is  approximately  2  °C  warmer.  Information  beyond 
SST  is  required  to  determine  which  case  is  closer  to  actual  conditions. 

Observations  of  vertical  temperature  profiles  provide  additional  insight  into  the  impact  of  ver¬ 
tical  coordinate  choice  on  subsurface  temperature.  Temperature  profile  data  are  obtained  from 
the  Marine  Environmental  Data  Service  (MEDS)  at  three  locations  available  during  1  January 
2000-31  March  2000.  MEDS  acquires,  processes,  quality  controls  and  archives  real-time  drifting 
buoy  messages  reporting  over  the  Global  Telecommunications  System  (GTS)  as  well  as  delayed 
mode  data  acquired  from  other  sources.  Further  information  about  the  MEDS  data  sets  is  avail¬ 
able  at  www.meds-sdmm.dfo-mpo.gc.ca/meds/Databases/Data_e.htm. 

Comparisons  of  subsurface  temperatures  from  NCOM  and  MEDS  are  made  at  three  locations 
(Fig.  12).  One  (43.9°N,  63.7°W)  is  in  the  western  North  Atlantic  near  the  shelf  break,  «200m 
depth  off  the  coast  of  Nova  Scotia.  The  second  (40.6°N,  130.1°W)  is  in  the  Japan/East  Sea 
(JES),  not  far  from  the  coast  but  in  fairly  deep  water  («2300  m),  and  the  third  is  in  the  deep,  open 
ocean  of  the  central  North  Pacific  (31.6°N,  173.20°W,  ~5000  m  depth).  These  samples  reflect  shelf 
edge,  semi-enclosed  sea  and  open  ocean  regions.  Temperature  profiles  at  these  locations  are  com¬ 
pared  on  5  February  2000,  15  March  2000  and  16  February  2000,  respectively.  Subsurface  tem¬ 
peratures  from  the  NCOM  models  are  extracted  for  the  same  dates  and  locations.  A  linear 
interpolation  is  applied  to  NCOM  data  to  match  depths  levels  of  the  MEDS  profiles,  allowing 
profiles  to  be  compared  at  the  same  depths. 

For  these  cases,  the  a-z  solutions  are  preferable  to  those  from  the  pure  z  simulation.  Solutions 
are  generally  similar  below  200  m,  a  zone  where  the  coordinate  systems  match.  Above  200  m,  the 
a-z  results  tend  to  be  in  closer  agreement  with  the  observations.  An  exception  is  the  JES 
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Fig.  12.  Comparisons  of  temperature  profiles  between  l/8°NCOM  and  MEDS  data  at  three  locations:  (left)  Shelf  edge 
off  Nova  Scotia,  43.9°N,  63.7°W,  200  m  depth,  5  February  2000;  (center)  Japan/East  Sea,  40.6°N,  130. 1°W,  2300  m 
depth,  15  March  2000;  (right)  central  North  Pacific,  31.6°N,  173.20°W,  5000  m  depth,  16  February  2000. 


observation,  where  the  z  simulation  has  smaller  temperature  errors  between  100  and  200  m.  How¬ 
ever  the  shape  and  associated  temperature  gradients,  more  important  than  actual  temperature  for 
acoustic  applications,  are  generally  better  from  the  a-z  simulation.  Neither  model  reproduces 
some  of  the  finer  scale  details,  such  as  the  warm  layer  in  the  upper  25  m  in  the  JES  sample  and 
the  temperature  inversion  at  a  depth  of  70  m  in  the  profile  near  Nova  Scotia.  Details  at  this  scale 
would  likely  require  assimilative  runs  at  higher  resolution  for  more  realistic  representation.  Dif¬ 
ferences  are  small  in  the  open  ocean  case,  as  expected. 


6,  Discussion  and  conclusions 

We  have  described  the  configuration  and  global  implementation  of  1/8°  global  NCOM.  The 
model  is  baroclinic,  hydrostatic,  Boussinesq,  and  has  a  free-surface.  It  uses  curvilinear  coordinates 
in  the  horizontal,  configured  with  a  rotated  reprojected  bipolar  grid  in  the  Arctic.  One  of  its  dis¬ 
tinguishing  characteristics  is  the  a-z  coordinate  in  the  vertical,  combining  z-levels  for  depths 
above  a  user-specified  depth,  here  137  m,  with  z  levels  below.  The  use  of  combined  a  and  z-level 
coordinates  in  a  single  model  provides  some  flexibility  in  setting  up  the  vertical  grid  to  address 
problems  that  are  expected  using  more  traditional  vertical  coordinates. 

A  pair  of  atmospherically-forced  free-running  simulations  conducted  with  the  model  using 
either  a-z  or  z-level  vertical  coordinates  demonstrate  the  flexibility  of  the  model  and  illustrate 
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some  consequences  arising  from  the  coordinate  choices.  It  is  noted  that  in  NCOM  z-level 
coordinates  the  bathymetry  is  rounded  to  the  nearest  model  level,  so  simulations  have  a  sub-opti¬ 
mal  representation  of  changes  in  the  bottom  depth.  The  model  solution  obtained  is  for  the  stair¬ 
step  bottom  used  in  the  model,  rather  than  the  actual  bathymetry.  Onshore  and  offshore 
barotropic  flows  can  be  noticeably  distorted  by  the  stairsteps  of  a  z-level  grid.  While  a  coordinates 
accurately  represent  the  changing  bottom  depth,  they  can  suffer  from  problems  with  their  horizon¬ 
tal  advection,  diffusion,  and  baroclinic  pressure  gradient  terms  in  regions  of  steep  bottom  slopes. 
Large  changes  in  depth  between  adjacent  grid  points  can  result  in  large  changes  in  the  model  fields 
between  grid  points  in  the  near-bottom  a  layers  due  to  large  vertical  gradients  in  the  oceanic  fields. 
This  can  present  problems  for  numerical  advection  schemes  that  require  gradients  to  be  resolved 
by  a  number  of  points  for  accuracy.  Horizontal  diffusion  between  adjacent  grid  points  that  are  at 
significantly  different  depths  can  effectively  result  in  vertical  diffusion  of  model  fields,  which  is  gen¬ 
erally  outside  of  turbulent  boundary  layers  where  oceanic  vertical  diffusion  tends  to  be  small. 

A  number  of  solutions  could  be  proposed  to  mitigate  shortcomings  of  various  vertical  coordi¬ 
nate  approaches.  One  solution  to  these  types  of  truncation  error  problems  is  to  increase  grid  res¬ 
olution.  However,  due  to  limited  computing  resources,  this  is  not  always  practical.  The  impact  of 
the  cost  of  more  expensive  numerical  methods  is  particularly  felt  in  global  scale  operational  mod¬ 
els  limited  by  operational  time  and  resource  constraints.  A  detailed  examination  of  a  wide  range 
of  mitigation  techniques  and  coordinate  definitions  is  beyond  the  scope  of  this  paper.  We  merely 
present  the  approach  used  in  global  NCOM,  the  hybrid  a-z  coordinate,  and  compare  samples  of 
its  results  with  a  parallel  case  using  a  more-traditional  z  coordinate.  The  combined  a-z  grid  allows 
the  use  of  multiple  a  coordinates  near  the  surface  to  take  up  changes  in  surface  elevation.  This 
appears  to  be  more  robust  at  longer  time  steps  than  a  single  free  surface.  In  shallow  water,  where 
flow  is  expected  to  be  more  sensitive  to  accurate  representation  of  bottom  depth,  the  a  coordinates 
allow  fidelity  to  bottom  topography  and  more  consistent  resolution  in  the  bottom  boundary  layer, 
z-levels  can  be  used  in  the  deeper  water,  e.g.,  below  the  continental  shelves  where  the  bottom 
slopes  can  be  very  steep,  to  avoid  some  of  the  problems  associated  with  a  coordinates  over  steep 
topography.  The  global  simulations  conducted  with  both  a-z  and  z-level  coordinates  to  assess  the 
extent  of  differences  caused  by  the  different  coordinate  systems  show  greatest  impact  in  nearshore 
regions  and  small  impact  in  open  ocean  regions,  as  expected.  For  global  NCOM,  comparisons 
with  temperature  observations  indicate  improved  overall  performance  using  the  a-z  coordinates 
relative  to  pure  z  coordinates. 

Operational  global  NCOM  provides  the  US  Navy  a  readily  available  first  look  at  mesoscale  cir¬ 
culation  conditions  in  any  region  of  the  global  ocean.  A  perhaps  more  important  objective  is  to 
provide  initial  and  boundary  conditions  for  fixed  and  relocatable  ocean  models.  Some  regional 
applications  have  been  tested  in  the  North  Atlantic,  East  Asian  Seas,  Intra-Americas  Seas,  Gulf 
of  Mexico,  and  US  west  coast  regions.  Analyses  of  these  applications  will  be  reported  in  future 
articles.  Further  objectives  of  global  NCOM  include  providing  capability  as  the  ocean  model 
component  for  a  global,  coupled,  air-ocean  modeling  system  and  hosting  an  embedded  Arctic 
ice  model.  NCOM  has  been  specifically  developed  to  be  suitable  for  coupling  in  a  regional  air/ 
ocean  system,  making  global  NCOM  a  candidate  for  use  as  the  ocean  component  in  a  global  cou¬ 
pled  system.  Work  on  embedding  the  Polar  Ice  Prediction  System  (PIPS)  3.0,  the  latest  generation 
of  the  Navy  ice  model,  is  scheduled  to  begin  in  late  2004.  Work  continues  on  refinements  to  the 
model  implementation  and  data  assimilation.  Additional  results  and  information  may  be  found  at 
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the  global  NCOM  web  page,  www.ocean.nrlssc.navy.mil/global_ncom.  An  evaluation  of  the  base¬ 
line  model  capability  is  necessary  to  complement  specialized  assessment  of  specific  applications. 
Capabilities  of  NCOM  over  the  global  ocean  are  examined  in  Kara  et  al.  (this  issue),  in  detail. 
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Appendix  A.  A  list  of  symbols 

Ah  horizontal  mixing  coefficient  for  temperature  or  salinity  (m2  s-1) 
A h0  constant  background  value  for  Ah  (10  m2  s-1) 

Am  horizontal  mixing  coefficient  for  momentum  (m2  s-1) 

^M0  constant  background  value  for  AM  (10  m2  s-1) 
cb  bottom  drag  coefficient 

/  Coriolis  parameter  (s-1) 

g  gravity  acceleration  (9.8  m  s-2) 

H  the  bottom  depth  (m) 

i  unit  vector  in  zonal  direction 

j  unit  vector  in  meridional  direction 

k  unit  vector  in  vertical  direction 

Kh  the  vertical  eddy  coefficient  for  temperature  or  salinity  (m2  s-1) 

Km  the  vertical  eddy  coefficient  for  momentum  (m2  s-1) 

p  pressure  (mb) 

<2r  net  solar  radiation  (W  m-2) 

S'  salinity  (psu) 

u  zonal  velocity  component  (ms-1) 

u  depth-averaged  zonal  velocity  (ms-1) 

v  meridional  velocity  component  (ms-1) 

v  velocity  vector 
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v  depth-averaged  meridional  velocity  (ms-1) 
t  time  (s) 

T  potential  temperature  (°C) 

x  distance  in  zonal  direction  (m) 

w  vertical  velocity  component  (ms-1) 

y  distance  in  meridional  direction  (m) 

z  distance  in  vertical  direction  (m) 

zb  distance  to  the  bottom  of  the  turbulent  layer  (m) 

zs  depth  at  which  the  grid  changes  from  a  to  z-level  (m) 

zt  distance  to  the  top  of  the  turbulent  layer  (m) 
y  a  function  describing  the  solar  extinction 
k  von  Karman  constant  (0.4) 

V  Laplacian  operator 

Vh  horizontal  Laplacian  operator 

p  water  density  (kg  m-3) 

p0  reference  water  density  (1025  kg  m-3) 
a  coordinate  transformation  (0  to  —1) 

6  latitude  (°) 

£  free  surface  elevation  above  the  undisturbed  value  at  z  =  0  (m) 


Appendix  B.  TAO  buoy  locations 

A  list  of  TAO  and  PIRATA  buoys  which  contain  daily-averaged  SST  data  used  for  NCOM 
SST  intercomparisons  in  2000.  Water  depths  at  buoy  locations  are  also  given.  It  should  be  noted 
that  daily  averages  are  constructed  from  hourly  SST  values.  Since  each  mooring  moves  in  time 
and  space  from  its  deployment  position,  we  calculated  average  position  based  on  the  historical 
latitude  and  longitude  data  for  each  buoy.  For  ease  of  notation,  nearest  integer  values  of  average 
latitude  and  longitude  is  used  for  each  buoy. 


TAO 

Location 

(Lat,  Lon) 

Depth  (m) 

Pacific 

(0°N,  095°W) 

(0.01  °N,  094.98°W) 

3318 

Pacific 

(0°N,  110°W) 

(0.05°N,  109.94°W) 

3935 

Pacific 

(0°N,  125°W) 

(0.20°S,  124.36°W) 

4780 

Pacific 

(0°N,  140°W) 

(0.05°N,  139.88°W) 

4359 

Pacific 

(0°N,  155°W) 

(0.05°N,  154.98°W) 

4647 

Pacific 

(0°N,  170°W) 

(0.04°S,  170.02°W) 

5616 

Pacific 

(0°N,  180°W) 

(0.02°N,  179.90°W) 

5439 

North 

(2°N,  095°W) 

(1.98°N,  095.03°W) 

3113 

North 

(2°N,  125°W) 

(1.97°N,  125.07°W) 

4745 

North 

(2°N,  140°W) 

(2.01  °N,  140.00°W) 

4387 

(continued  on  next  page ) 
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TAO 

Location 

(Lat,  Lon) 

Depth  (m) 

North 

(2°N,  155°W) 

(2.00°N,  154.97°W) 

4661 

North 

(2°N,  156°E) 

(2.00°N,  156.02°E) 

2594 

North 

(2°N,  170°W) 

(2.00°N,  170.03°W) 

5408 

North 

(2°N,  180°W) 

(2.02°N,  179.80°W) 

5477 

South 

(2°S,  110°W) 

(2.00°S,  109.98°W) 

3950 

South 

(2°S,  125°W) 

(2.04°S,  124.90°W) 

4751 

South 

(2°S,  140°W) 

(2.02°S,  139.95°W) 

4320 

South 

(2°S,  155°W) 

(1.99°S,  154.96°W) 

4988 

South 

(2°S,  170°W) 

(2.17°S,  170.04°W) 

4983 

South 

(2°S,  1 80°W) 

(2.01  °S,  179.90°W) 

5353 

North 

(5°N,  110°W) 

(4.94°N,  109.99°W) 

3991 

North 

(5°N,  125°W) 

(5.32°N,  124.94°W) 

4395 

North 

(5°N,  155°W) 

(4.95°N,  155.09°W) 

4597 

North 

(5°N,  156°E) 

(4.99°N,  156.06°E) 

3607 

North 

(5°N,  165°E) 

(5.05°N,  164.98°E) 

4780 

North 

(5°N,  180°W) 

(4.99°N,  179.93°W) 

5680 

South 

(5°S,  095°W) 

(5.02°S,  095.07°W) 

3831 

South 

(5°S,  110°W) 

(5.00°S,  109.99°W) 

3605 

South 

(5°S,  125°W) 

(5.02°S,  124.95°W) 

4561 

South 

(5°S,  140°W) 

(5.01  °S,  139.90°W) 

4359 

South 

(5°S,  155°W) 

(4.99°S,  154.98°W) 

5014 

South 

(5°S,  165°E) 

(5.00°S,  165.20°E) 

5418 

South 

(5°S,  180°W) 

(4.98°S,  179.90°W) 

5664 

North 

(8°N,  155°W) 

(7.97°N,  1 55.00°W) 

5249 

North 

(8°N,  156°E) 

(7.97°N,  156.03°E) 

4920 

North 

(8°N,  165°E) 

(7.98°N,  165.06°E) 

5210 

North 

(8°N,  170°W) 

(8.01  °N,  170.02°W) 

5553 

South 

(8°S,  110°W) 

(8.05°S,  109.93°W) 

3465 

South 

(8°S,  125°W) 

(7.96°S,  125.02°W) 

4564 

South 

(8°S,  155°W) 

(8.26°S,  155.02°W) 

5341 

Atlantic 

(0°N,  023°W) 

(0.02°N,  23.12°W) 

3843 

Atlantic 

(0°N,  035°W) 

(0.37°S,  35.06°W) 

4421 

North 

(15°N,  038°W) 

(14.94°N,  38.39°W) 

5382 

Appendix  C.  NODC  buoy  locations 

The  same  as  Appendix  B  but  for  the  NODC  buoys.  The  5-digit  numbers  column  represents 
buoy  ID  numbers  used  by  NODC. 
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NODC 

Location 

(Lat,  Lon) 

Depth  (m) 

(17N,  153°W) 

(17.43N,  152.52°W) 

5303 

(19N,  161°W) 

(19.17N,  160.72°W) 

4943 

(23N,  162°W) 

(23.40N,  162.25°W) 

3257 

(26N,  090°W) 

(25.92N,  089.67°W) 

3246 

(26N,  094°W) 

(25.88N,  093.57°W) 

3200 

42041 

(27N,  090°W) 

(27.22N,  090.43°W) 

1435 

(29N,  079°W) 

(28.88N,  078. 52° W) 

841 

(29N,  080°W) 

(28.50N,  080.18°W) 

42 

(29N,  085°W) 

(28.50N,  084.50°W) 

53 

(29N,  086°W) 

(28.78N,  086.03°W) 

283 

(29N,  088°W) 

(29.20N,  088.20°W) 

237 

(29N,  094°W) 

(29.23N,  094.40°W) 

15 

(30N,  089°W) 

(29.55N,  088.50°W) 

40 

(3 IN,  081°W) 

(31.40N,  080.87°W) 

18 

(32N,  120°W) 

(32.43N,  119.52°W) 

1393 

(33N,  079°W) 

(32.50N,  079.10°W) 

36 

(34N,  119°W) 

(33.73N,  119.08°W) 

859 

(34N,  120°W) 

(34.23N,  119.83°W) 

417 

(34N,  121°W) 

(34.25N,  120.65°W) 

598 

(35N,  121°W) 

(34.70N,  120.97°W) 

384 

(37N,  122°W) 

(36.75N,  122.42°W) 

1920 

(38N,  075°W) 

(38.45N,  074.70°W) 

28 

(38N,  123°W) 

(38.22N,  123.32°W) 

122 

(38N,  130°W) 

(37.98N,  129.98°W) 

4599 

(39N,  124°W) 

(39.22N,  123.95°W) 

264 

(40N,  073°W) 

(40.25N,  073.17°W) 

40 

(40N,  125°W) 

(40.42N,  124.52°W) 

82 

(41N,  069°W) 

(40.50N,  069.42°W) 

62 

(42N,  071°W) 

(42.35N,  070.68°W) 

55 

(42N,  124°W) 

(41.85N,  124.37°W) 

47 

(43N,  069°W) 

(42.88N,  068.93°W) 

21 

(44N,  070°W) 

(43.52N,  070.13°W) 

18 

(45N,  125°W) 

(44.62N,  124.52°W) 

130 

(46N,  131°W) 

(46.05N,  131.02°W) 

2779 

(47N,  125°W) 

(47.33N,  124.75°W) 

132 

(56N,  148°W) 

(56.28N,  148.17°W) 

4206 

46035 

(57N,  178°W) 

(56.90N,  177.80°W) 

3662 

46061 

(60N,  147°W) 
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